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Abstract 

We propose new measures of localization and cooperativity for the anal- 
ysis of atomic rearrangements. We show that for both clusters and bulk 
material cooperative rearrangements usually have significantly lower bar- 
riers than uncooperative ones, irrespective of the degree of localization. 
We also find that previous methods used to sample stationary points are 
biased towards rearrangements of particular types. Linear interpolation 
between local minima in double-ended transition state searches tends to 
produce cooperative rearrangements, while random perturbations of all the 
coordinates, as sometimes used in single-ended searches, has the opposite 
effect. 
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1 Introduction 



The potential energy surface (PES) governs the observed structure, dynamics 
and thermodynamics of any molecular system. It is often possible to gain new 
insight into these properties by expressing them in terms of stationary points 
of the PES, i.e. points where the gradient of the potential vanishes. 1 The most 
important stationary points are minima and the transition states that connect 
them. Here we define a minimum as a stationary point where the Hessian, the 
second derivative matrix, has no negative eigenvalues, while a transition state is 
a stationary point with precisely one such eigenvalue. 2 

The number of stationary points on the PES generally scales exponentially 
with system size, 3-7 which necessitates an appropriate sampling strategy of some 
sort for larger systems. In particular, to analyse dynamical properties a database 
of local minima and the transition states that connect them is usually constructed, 
which generally involves extensive use of single-ended and double-ended transi- 
tion state searching techniques (see Refs. 1,8,9 and references therein). Single- 
ended transition state searches only require an initial starting geometry. However, 
double-ended searches require two endpoint geometries, a mechanism to gener- 
ate a set of configurations between them, and a suitable functional (or gradient) 
to be evaluated and minimised. The most successful single- and double-ended 
methods currently appear to be based upon hybrid eigenvector-following 10-22 and 
the nudged elastic band approach, 8, 9,23-27 respectively. The two search types are 
often used together, since double-ended transition state searches do not produce 
a tightly converged transition state and further refinement may be needed. 1,8 

Any path connecting two minima on PES can be broken down into elemen- 
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tary rearrangements, each of which involves a single transition state. The cor- 
responding mechanism can be analyzed in detail by calculating the two unique 
steep est- descent paths that lead downhill from the transition state. 

The number of elementary rearrangements, as defined above, increases expo- 
nentially with system size as for the number of transition states. For instance, 
there are approximately 30,000 such pathways on the PES of the 13-atom cluster 
bound by the Lennard- Jones potential. When permutation-inversion isomers are 
included, this number increases by a factor of order 2 x iV!. 1 Two activation barri- 
ers can be defined for each pathway in terms of the energy difference between the 
transition state and each of the minima. For non-degenerate rearrangements 1,28 
the two sides of the path are termed uphill and downhill, where the uphill bar- 
rier is the larger one, which leads to the higher minimum. The barriers and the 
normal modes of the minima and transition states can be used to calculate rate 
constants using harmonic transition state theory. 29 " 31 

For each local minimum a catchment basin can be defined in terms of all the 
configurations from which steepest-descent paths lead to that minimum. 32 Some 
of these paths originate from transition states on the boundary of the catchment 
basin, which connect a given minimum to adjacent minima. The integrated path 
length for such rearrangements provides a measure of the separation between 
local minima, and may be related to the density of stationary points in configu- 
ration space. The integrated path length is usually approximated as the sum of 
Euclidean distances between configurations sampled along appropriate steepest- 
descent paths. 1 It provides a convenient coordinate for monitoring the progress 
of the reaction. 

Calculated pathways can always be further classified mechanistically. For ex- 
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ample, some rearrangements preserve the nearest-neighbour coordination shell for 
all the atoms. In previous studies of bulk models these cage-preserving pathways 
were generally found to outnumber the more localized cage-breaking processes, 
which are necessary for atomic transport. 33 It was found that the barriers for 
cage-breaking and cage-preserving processes were similar for bulk LJ systems, 
while the cage-breaking mechanisms have significantly higher barriers for bulk 
silicon modelled by the Stillinger- Weber potential. 33 

For minima separated by increasing distances in configuration space, the path- 
ways that connect them are likely to involve more and more elementary steps, 
and are not unique. Finding such paths in high-dimensional systems can become 
a challenging task. 8 ' 34 Some difficulties have been attributed to instabilities and 
inefficiencies in transition state searching algorithms, 8,35 as well as the existence 
of very different barrier and path length scales. 8 A new algorithm for locating 
multi-step pathways in such cases has recently been proposed. 8,34 

In the present work we have used a doubly-nudged elastic band (DNEB) 
method 8 in conjunction with eigenvector-following (EF) algorithms 10-22 to locate 
rearrangement pathways in various systems. The Lennard- Jones (LJ) potential, 
defined as V(r) = 4e [(a/r) 12 — (er/r) 6 ], where r is the distance between two 
atoms, e is the depth of the potential energy well, and 2 1 / 6 <r is the pair equilibrium 
separation, was used to describe the 13- and 75-atom Lennard- Jones clusters, LJ13 
and LJ75. We have also considered a binary LJ (BLJ) system with parameters 
a AA = 1.0, a BB = 0.88, a AB = 0.8, e AA = 1.0, e BB = 0.5, e AB = 1.5, where 
A and B are atom types. The mixture with A : B ratio 80 : 20 provides a 
popular model bulk glass-former. 33,36 We employed a periodically repeated cubic 
cell containing 205 A atoms and 51 B atoms. The density was fixed at 1.2a^ A 
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and the Stoddard- Ford scheme was used to prevent discontinuities. 3 

The motivation for this paper was our observation that construction of some 
multi-step pathways using the connection algorithm described in Ref. 8 is partic- 
ularly difficult. We were unable to relate these difficulties to simple properties of 
the pathways such as the integrated path length, the uphill and downhill barriers, 
or the barrier and path length asymmetries. Instead, more precise measures of 
localization and cooperativity are required, as shown in the following sections. It 
also seems likely that such tools may prove useful in analysing the dynamics of 
supercooled liquids, where processes such as intrabasin oscillations and interbasin 
hopping have been associated with different rearrangement mechanisms. 38 In par- 
ticular, cooperativity is believed to play an important role at low temperatures 
in glass-forming systems, 39 and dynamic heterogeneity may result in decoupling 
between structural relaxation and transport properties for supercooled liquids. 40 

2 Methods 
2.1 Localization 

The outcome of a pathway calculation for an atomic system will generally be a set 
of intermediate geometries, and the corresponding energies, for points along the 
two unique steepest-descent paths that link a transition state to two local minima. 
This discrete representation is a convenient starting point for our analysis of 
localization and cooperativity. We number the structures along the path j = 
1, 2, . . . , Nf starting from one of the two minima and reversing the other steepest- 
descent path, so that structure Nf corresponds to the other minimum. The 
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transition state then lies somewhere between frames 1 and Nf. We define the 
three-dimensional vector Ti(j) to contain the Cartesian coordinates of atom % for 
structure j, i.e. r = ^Xj(j), Yi(j), Zi(j) J , where Xi(J) is the X coordinate of 
atom % in structure j, etc. 

For each atom % we also define the displacement between structures j — 1 and 
j as 

diU) = \nU)-*iU-i)[ (!) 

Hence the sum of displacements 

<*i = X>0'), (2) 

i=2 

is an approximation to the integrated path length for atom i, which becomes 
increasingly accurate for smaller step sizes. The total integrated path length, s, 
is approximated as 



N 



N f 

S 

i=2 

where N is the total number of atoms, s is a characteristic property of the 



,5>^') 2 ' ( 3 ) 
\ i=i 



complete path, and is expected to correlate with parameters such as the curvature 
and barrier height for short paths. 1 ' 41 ' 42 

The set {d\, ofo, ■ ■ ■ , c?v} containing all N values of di will be denoted {d}, and 
analogous notation will be used for other sets below. We will also refer to the 
frequency distribution function, which can provide an alternative representation 
of such data. 43 For example, the frequency distribution function F for a given 
continuous variable, x, tells us that x occurs in a certain interval F(x) times. 

Our objective in the present analysis is to provide a more detailed descrip- 
tion of the degree of 'localization' and 'cooperativity' corresponding to a given 
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pathway. The first index we consider is N p , which is designed to provide an es- 
timate of how many atoms participate in the rearrangement. We will refer to 
a rearrangement as localized if a small fraction of the atoms participate in the 
rearrangement, and as delocalized in the opposite limit. The second index we 
define, N c , is intended to characterize the number of atoms that move simultane- 
ously, i.e. cooperatively. We will refer to a rearrangement as cooperative if most 
of the atoms that participate in the rearrangement move simultaneously, and as 
uncooperative otherwise. 

The nth moment about the mean for a data set {x\,X2, ...,xm} is the expec- 
tation value of (xi — (x)) n , where (x) = YltLi x i/M and M is the number of 
elements in the set. Hence for the set {d} defined above we define the moments, 
m n , as m n = Yli=i — ( c 0)" /N. The kurtosis of the set {d} is then defined as 
the dimensionless ratio 7(d) = m 4 /(m 2 ) 2 , and provides a measure of the shape 
of the frequency distribution function corresponding to {d}. If only one of the 
atoms moves, or all atoms except one move by the same amount, then 7(d) = N. 
Alternatively, if half the atoms move by the same amount whilst the others are 
stationary, then 7(d) = 1. Hence, a distribution with a broad peak and rapidly 
decaying tails will have a small kurtosis, 7 ~ 0(1), while a distribution with 
a sharp peak and slowly decaying tails will have a larger value (Fig. ^). The 
kurtosis can therefore identify distributions that contain large deviations from 
the average value. 43 For comparison, a Gaussian distribution has 7 = 3 and a 
uniform distribution has 7 = 1.8. 

The above results show that the kurtosis 7(d) can be used to quantify the 
degree of localization or derealization of a given rearrangement. However, it has 
the serious disadvantage that highly localized and delocalized mechanisms both 
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have large values of 7(d). Since we are interested in estimating the number of 
atoms that move relative to the number with small or zero displacements, a better 
approach is to use moments taken about the origin, rather than about the mean, 
i.e. m' n = ^idy/N, following Stillinger and Weber. 44 Note that while mi = 0, 
the first moment m! Y is the mean value. We therefore estimate the number of 
atoms that participate in the rearrangement, N p , as 



where ^y'(d) = m'J (m' 2 ) 2 . For the system with N atoms, if only one atom moves 
N p — 1, while if K atoms move by the same amount, N p = K. 

A similar index to N p has been employed in previous work 1 ' 17,45 using only the 
displacements between the two local minima, which corresponds to taking Nf = 2 
in Eq. (|2j). Using di values based upon a sum of displacements that approximates 
the integrated path length for atom i, rather than the overall displacement be- 
tween the two minima, better reflects the character of the rearrangement, as it 
can account for the nonlinearity of the pathway. To describe this property more 
precisely we introduce a pathway nonlinearity index defined by a = (s — D)/s, 
where D is the Euclidean distance between the endpoints, 



We calculated the a values for a database of 31,342 single transition state 
pathways of LJ 75 (hereafter referred to as the LJ 75 database). The average value 
of a was 0.4 with a standard deviation of 0.2, and, hence, there is a significant 
loss of information if 7' is calculated only from the endpoints using Nf = 2. 
Comparison of the two indices for the LJ 75 database revealed many examples 
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where neglect of intermediate structures produces a misleading impression of the 
number of atoms that move. The definition in equation (@J) is therefore suggested 
as an improvement on previous indices of localization. 1 ' 17,44 ' 45 

2.2 Cooperativity 

N p atoms can participate in a rearrangement according to a continuous range 
of cooperativity. At one end of the scale there are rearrangements where N p 
atoms all move simultaneously [see Fig. Efa)]- Although these paths exhibit the 
highest degree of correlated atomic motion they do not usually pose a problem for 
double-ended transition state search algorithms. 8,9 Linear interpolation between 
the two minima tends to generate initial guesses that lie close to the true pathway, 
particularly if a ~ 0. At the opposite extreme, atoms can move almost one at a 
time, following a 'domino' pattern [see Fig. E£b)]. Locating a transition state for 
such rearrangements may require a better initial guess, since linear interpolation 
effectively assumes that all the coordinates change at the same rate. 

The degree of correlation in the atomic displacements can be quantified by 
considering the displacement 'overlap' 



where the index k indicates that O was calculated for k atoms numbered c(l), 
c(2),...,c(k). c is a k- dimensional vector that represents a particular choice of k 
atoms from N, and hence there are C% = N\/k\(N — k)\ possible values of Ok- 
The index O can be thought of as a measure of how the displacements of the 
atoms c(l), c(2), etc. overlap along the pathway. For example, if two atoms move 
at different times then O2 is small for this pair because the minimum displacement 
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in Eq. (jUJ) is always small. However, if both atoms move in the same region of 
the path then O2 is larger. 

We now explain how the statistics of the overlaps, Ok, can be used to extract 
a measure of cooperativity (Fig. EJ). Suppose that m atoms move simultaneously 
in a hypothetical rearrangement. Then all the overlaps Ok for k > m will be 
relatively small, because one or more atoms are included in the calculation whose 
motion is uncorrelated with the others. For overlaps Ok with k ^ m the set of 
Ok for all possible choices of k atoms from N will exhibit some large values and 
some small. The large values occur when all the chosen atoms are members of the 
set that move cooperatively, while other choices give small values of Ok- Hence 
the kurtosis of the set {Ok}, '-/'(Ok), calculated from moments taken about the 
origin, will be large for k ^ m, and small for k > m. To obtain a measure of 
how many atoms move cooperatively we could therefore calculate ^'(02), l'{0^), 
etc. and look for the value of k where 7'(Ofc) falls in magnitude. However, to 
avoid an arbitrary cut-off, it is better to calculate the kurtosis of the set {^'(02), 
^'{Os), 7'(Ofc)}, or 7'[7'(0)] for short. There are N — 2 members of this set, and 
by analogy with the definition of N p = N/j'(d), we could define a cooperativity 
index N c = (N - 2)/i[i(0)\ + 1. Then, if i{0 2 ) is large, and all the other 
J 1 (Ok) are small, we obtain J , ('Y(0)) ~ N — 2 and N c ~ 2, correctly reflecting 
the number of atoms that move together. 

In practice, there are several problems with the above definition of iV c . Cal- 
culating iV c in this way quickly becomes costly as the number of atoms and/or 
number of frames in the pathway increases, because the number of elements in 
the set {Ok} varies combinatorially with k. Secondly, as k approaches iV the dis- 
tribution of all the possible values for Ok becomes more and more uniform. Under 
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these circumstances deviations from the mean that are negligible in comparison 
with the overall displacement can produce large kurtosises. Instead, we suggest 
a modified (and simpler) definition of N c , which better satisfies our objectives. 

We first define the overlap of atomic displacements in a different manner. It 
can be seen from equation (JHJ) that the simultaneous displacement of / atoms is in- 
cluded in each set of overlaps {Ok} with k ^ I. For example, if three atoms move 
cooperatively then both the {O2} and {O3} sets will include large elements cor- 
responding to these contributions. Another redundancy is present within {Ok}, 
since values in this set are calculated for all possible subsets of k atoms and the 
displacement of each atom is therefore considered more than once. However, we 
can avoid this redundancy by defining a single k— overlap, rather than dealing 
with C% different values. 

Recall that di(j) is the displacement of atom % between frames j — 1 and j. The 
ordering of the atoms is arbitrary but remains the same for each frame number 
j. We now define Aj(j) as the displacement of atom % in frame j, where index % 
numbers the atoms in frame j in descending order, according to the magnitude 
of di(j), e.g. atom 1 in frame 2 is now the atom with the maximum displacement 
between frames 1 and 2, atom 2 has the second largest displacement etc. As the 
ordering may vary from frame to frame, the atoms labelled % in different frames 
can now be different. This relabelling greatly simplifies the notation we are about 
to introduce. Consider the /c-overlap defined as 



where k ranges from 1 to N, A to t = Ylj=2^U) anc ^ ^N+i(j) is defined to be 
zero for all j. A schematic illustration of this construct is presented in Fig.0] For 
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example, if only two atoms move in the course of the rearrangement, and both 
are displaced by the same amount (which may vary from frame to frame), the 
only non-zero overlap will be 02- 

We can now define an index to quantify the number of atoms that move 
cooperatively as 

N 

Nc = Y^kQ k . (8) 

fe=i 

If only one atom moves during the rearrangement then N c = 1, while if K atoms 
displace cooperatively during the rearrangement then N c = K. This definition 
is independent of the total displacement, the integrated path length, and the 
number of atoms, which makes it possible to compare N c indices calculated for 
different systems. 

3 Results 

Fig.Elshows results for the most cooperative and uncooperative processes we have 
found for the LJ 75 cluster that are localized mainly on two atoms. In these calcu- 
lations we have used the database of transition states that was found previously as 
the result of a discrete path sampling calculation conducted for this system. 34 ' 46 
The cooperative rearrangement [Fig. |2fa,c)] is the one with the maximum two- 
overlap 02- For this pathway 2 = 0.7, N p = 3.4, and N c = 7.7. The values of 
N p and iV c both reflect the fact that the motion of the two atoms is accompanied 
by a slight distortion of the cluster core. This example shows that while N p and 
N c allow us to quantify localization and cooperativity, and correctly reflect the 
number of atoms that participate and move cooperatively in ideal cases, there will 
not generally be a simple correspondence between their values and the number of 
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atoms that move. This complication is due to the fact that small displacements of 
atoms in the core will generally occur, no matter how localized the rearrangement 
is. In addition, the data reduction performed in equations (J1J and (jBJ) means that 
a range of pathways can give the same value for N p or N c . Since the size of the 
contribution from a large number of small displacements depends on the shape 
of the displacement distribution function the number of possibilities grows with 
the size of the index. 

The uncooperative rearrangement depicted in the Fig. El^b,d) was harder to 
identify. In principle we could have selected the pathway that maximizes 0i from 
all the rearrangements localized on two atoms. However, this approach picks out 
rearrangements localized on one atom, where distortion of the core accounts for 
the value of N p > 1. Instead, we first selected from all the rearrangements 
with N p < 4 those where two atoms move by approximately the same amount, 
while the displacement of any other atom is significantly smaller. These are 
the rearrangements that maximize l/7({di, e^}), where d\ and d>2 are the total 
displacements of the two atoms that move the most. After this procedure we 
selected the rearrangement with the maximum value of Q±. Fig. Efb) shows that 
this rearrangement features the displacement of one atom at a time, and the 
atom that moves first also moves last. For this pathway the values of Gi, N p 
and N c are 0.7, 3.8 and 5.3, respectively. Further illustrations and movies of the 
corresponding rearrangements are available online. 47 

Fig. El illustrates several general trends that we have observed for cluster re- 
arrangements. Firstly, we have found that the barrier height is smaller for the 
cooperative rearrangements [Fig. Efc,d)]. Usually atoms that move cooperatively 
are neighbours. Rearrangements generally involve a change of the environment 
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for the atoms that move. Cooperative motion can reduce this perturbation since 
for any of the participants the local environment is partly preserved because it 
moves with the atom in question. Flatter points on the energy profile [circled in 
Fig.|2fd)] usually signify a change in the mechanism, i.e. one group of atoms stops 
moving and another group starts. By comparing (b) and (d) in Fig.|2]we conclude 
that flatter points on the energy profile correlate with the most cooperative parts 
of this rearrangement. 

A simple correlation between barrier heights and N p and iV c does not seem 
to exist. The barrier height is not a function of cooperativity alone, but also 
of the energetics of the participating atoms. The way the N p and N c indices 
have been defined can make them insensitive to details of the rearrangements 
that will affect the energetics. For instance, neither index depends on the lo- 
cation of the participating atoms or the directionality of their motion. In most 
cases cooperatively moving particles are adjacent, i.e. localized in space; however, 
long-distance correlations of atomic displacements also occur. One such case is 
depicted in Fig. |3] (a). This path is nearly symmetric with respect to the inte- 
grated path length (it = 0.01, see Ref. 8), but is very asymmetric with respect to 
the uphill and downhill barrier heights ((3 = 0.91, see Ref. 8). This rearrangement 
has N p = N c = 10. Interestingly, N p and N c calculated separately for both sides 
of the pathway are very similar, i.e. the two steepest-descent pathways cannot be 
distinguished using these indices. Close inspection of this rearrangement reveals 
that one side of the pathway involves the rearrangement of two atomic triplets 
that share a vertex, while the other side involves the drift of all five atoms on 
the surface of the cluster [see the insets in Fig. (a)]. Although N c does not 
distinguish these cases, the motion in the second side of the path is more coop- 
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erative. The participating atoms move together, which results in a significantly 
lower downhill barrier. 

N p and N c also describe properties of the whole pathway. A significant num- 
ber of pathways that we observed were rather non-uniform, i.e. very cooperative 
phases alternated with uncooperative ones. To distinguish such pathways in the 
LJ75 database we calculated a set {N p } containing N p 'b evaluated for each pair of 
adjacent frames. Then a selection of pathways was made with m2/(m' 1 ) 2 < 0.01, 
where rri2/(m' 1 ) 2 is the moment ratio evaluated for the set {N p }. While this proce- 
dure ensured that N p corresponds closely to the number of atoms that moves be- 
tween any two snapshots of the rearrangement, it did not distinguish cases where 
different atoms contribute to the value of N p in different frames [see Fig. Efb)]. 
The average uphill and downhill barriers for this subset of rearrangements are 100 
times smaller than the average barriers for the complete LJ 75 database (Table Q]). 
Fig. El shows that N p and N c calculated for these rearrangements are highly cor- 
related and span a range of values, implying that widely different pathways are 
represented. Finally, all the selected pathways are an order of magnitude shorter 
than the average path length for the whole database, even though this database 
contains many short rearrangements localized on one atom. 

Fig. 0shows the values of the N p and iV c indices plotted against each other for 
pathway databases calculated for LJ13, LJ75 and BLJ256. The two databases for 
BLJ 256 labelled as 1 and 12 are taken from Ref. 33 and correspond to databases 
BLJ1 and BLJ12 in that paper. BLJ1 and BLJ12 were obtained using two dif- 
ferent sampling schemes intended to provide an overview of a wide range of con- 
figuration space and a thorough probe of a smaller region, respectively. These 
databases were constructed by systematic exploration of the PES, and we refer 
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the reader to the original work for further details. 33 Each BLJ 2 56 database con- 
tains 10, 000 transition states. The LJ 13 and LJ 75 databases were obtained in 
discrete path sampling (DPS) studies 34 ' 46 and contain 28, 756 and 31, 342 transi- 
tion states, respectively. Fig. [7| is a density plot where darker shading signifies a 
higher concentration of data points. The outlying points are connected by a solid 
line to define the area in which all the points lie. Fig. shows that as N p grows 
the allowed range of N c increases, especially for LJ 13 . For the LJ75 database rear- 
rangements with iV c > N/2 appear to be very rare or poorly sampled. Fig. also 
shows that for all these systems rearrangements localized on two or three atoms 
dominate. This result may be an intrinsic property. However, it may also be due 
to the geometric perturbation scheme used in producing the starting points for the 
transition state searches employed in generating these databases. For databases 
LJ 13 , BLJ1 and BLJ12 there are significantly more rearrangements with larger 
values of N p and N c compared to LJ75, which suggests that the abundance of 
very localized rearrangements for clusters may be a surface effect. 

Fig. |H] depicts the average barrier as a function of the participation and co- 
operativity indices. N p and N c were calculated separately for both sides of the 
pathway and the corresponding barriers (uphill or downhill) were averaged to 
produce a density plot of barrier height. This figure illustrates that for each 
system cooperative rearrangements have the lowest barriers, irrespective of the 
value of N p . For clusters, cooperative rearrangements have lower barriers than 
uncooperative rearrangements with N p as small as 1 — 3, while for bulk barriers 
corresponding to rearrangements with low N p become comparable to these for 
very cooperative rearrangements with high N p . 

In further computational experiments we find that attempts to connect the 
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endpoints of uncooperative pathways using the algorithm described in Ref. 8 
either required more images and iterations or converged to an alternative pathway. 
In some cases additional difficulties arose, such as convergence to a higher index 
saddle instead of a transition state, which can happen if the linear interpolation 
guess conserves a symmetry plane. Fig. |H1 shows N c calculated from Eq. (JHJ) 
plotted against N p for the LJ75 pathway database. Knowing the integrated path 
length, s, for each pathway we started doubly nudged elastic band calculations 
with three images per unit of distance and 30 iterations per image. Most of the 
points that correspond to runs that failed or converged to an alternative pathways 
are concentrated in the region of small values for N c /N p . 

As can be seen from Fig. Q the LJ13 database contains significantly more 
pathways with large values of N p and N c compared to LJ75, where most of the 
rearrangements are localized on two atoms. The fact that the LJ 13 database 
is almost exhaustive then suggests that localized rearrangements either start to 
dominate as the system size increases or that the sampling scheme used for LJ75 
was biased towards such mechanisms. Systematic sampling of the configuration 
space for stationary points often employs perturbations of every degree of free- 
dom followed by minimization. 48 The LJ13 database was obtained in this fashion, 
while the LJ 75 database was generated during the discrete path sampling (DPS) 
approach. 34 In this procedure discrete paths are perturbed by replacing local 
minima with structures obtained after perturbing all the coordinates and mini- 
mizing. To investigate whether the perturbation scheme can affect the resulting 
database of stationary points in more detail we consider the case of LJ13, since 
nearly all the transition states are known. Fig. EH presents the results of two 
independent runs aimed at locating most of the transition states for this system. 
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Every cycle a perturbation was applied to a randomly selected transition state 
from the database and the resulting geometry was used as a starting point for 
a new transition state search using eigenvector-following. 17 ' 20-22,49 Only distinct 
permutation-inversion isomers were saved. In the first run (bottom curve) every 
degree of freedom was perturbed by OAx, where x is a random number in the 
interval [—1,1]. 48 For the second run (top curve) we introduced a perturbation 
scheme including correlation. 2 «C K $C N/2 atoms out of N were displaced by 
a vector 0.4(xi, X2, £3), where the components X\, X2 and X3 are again random 
numbers drawn from [—1, 1]. The K atoms to be displaced were selected based 
on their relative positions in the cluster. One atom was first selected at random, 
while the remaining K — 1 were chosen to be its closest neighbours. The top 
curve was generated from a run with K = 6. Both runs required approximately 
the same time to produce two nearly identical databases, each containing about 
29, 000 pathways. However, as can be seen from Fig. random perturbation of 
all the degrees of freedom results in uncooperative rearrangements being found 
first, while employing correlated perturbations has the opposite effect. 

4 Conclusions 

The most important result of this work is probably the introduction of an index 
to quantify the cooperativity of atomic rearrangements. With this new measure 
it becomes possible to correlate cooperativity and barrier heights, and to show 
that cooperative rearrangements generally have lower barriers and shorter path 
lengths. We hope that these results will shed new light on relaxation mechanisms 
in complex systems, such as glasses and biomolecules, in future applications. For 
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example, in a peptide or protein a large geometrical change can result from a 
rearrangement that could be described in terms of a single dihedral angle. In 
glasses and supercooled liquids an important research goal is to understand how 
observed dynamical properties, such as atomic diffusion and correlation func- 
tions, 33 ' 38 ' 50 ' 51 are related to features of the underlying PES. The classification 
of elementary rearrangements as 'cage-breaking' or 'cage-preserving', 33,52 and the 
emergence of structures such as 'megabasins' 33, 52-54 can now be investigated more 
precisely in terms of localization and cooperativity. 

We have also demonstrated that cooperative rearrangements are relatively 
easy to characterize using double-ended transition state searching algorithms, 
since linear interpolation produces an effective initial guess. Uncooperative rear- 
rangements are usually harder to find using such methods, and alternative initial 
guesses may be helpful in these cases. 

Single-ended transition state searching has been used both in conjunction 
with double-ended methods, and as a way to sample potential energy surfaces for 
stationary points. Stationary point databases constructed using random pertur- 
bations followed by quenching are likely to be biased towards uncooperative rear- 
rangements. We have therefore outlined a strategy for generating initial guesses 
appropriate to single-ended transition state searching algorithms, which instead 
favours cooperative rearrangements. This approach also includes a parameter 
that is likely to influence the degree of localization. 
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6 Tables 



Table 1: Average uphill and downhill barriers, average integrated path length for 
the LJ75 rearrangement pathway database. Values are given for the whole database 
containing 31,342 paths, and for a subset containing the 57 most cooperative paths. 
The units of energy and distance are e and a, respectively. 

All Cooperative 
Uphill barrier 3.03 0.06 
Downhill barrier 0.97 0.03 
Path length 3.08 0.58 
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Figure Captions 

Two frequency distribution functions F\ and F2 of a continuous variable 
x are contrasted. Both functions have the same average m! 1 and standard 
deviation \frn 2 . However, due to the long tails, F\ has a significantly larger 
fourth moment 7714 and, hence, a larger kurtosis, 7. 

Comparison of cooperative (a) and uncooperative (b) rearrangements of 
the LJ 75 cluster, for mechanisms that are localized mainly on two atoms. 
The displacement d as a function of the frame number j is shown for the 
two atoms that move the most. Panels (c) and (d) illustrate the potential 
energy V as a function of frame number j for the rearrangements in (a) and 
(b), respectively. Dashed circles indicate flatter parts of the energy profile, 
which correspond to the most cooperative regions of the pathway. 

V(02) plotted against ^'(d) for the LJ 75 pathway database. The figure 
shows how ^'{O-z) can discriminate between rearrangements that have sim- 
ilar values of ^'(d) but different cooperativity. The data point for the most 
cooperative rearrangement localized on two atoms has N p = N c = 2. 

The G indices for a hypothetical rearrangement localized on three atoms. 
For each of these atoms the displacement d as a function of frame number 
j is shown. The values in successive frames are connected with dotted 
(di), dashed (g^) and solid (d^) lines. The corresponding contributions to 
0i, ©2, and 03 are shown as shaded squares and are labelled accordingly. 
If the remaining N — 3 atoms do not participate, and the area of one square 
is S, the only non-zero overlaps will be ©i, G2, and 63 with values 5/9, 3/9 
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and 1/9, respectively. 

5. Two limitations of the cooperativity index N c . (a) N c is not sensitive to the 
spatial positions of the cooperatively moving atoms, nor to the directionality 
of their motion. The energy profile is depicted for a rearrangement of the 
LJ75 cluster, which is very asymmetric with respect to barrier height but has 
similar integrated path lengths on either side of the transition state. The 
cooperativity index N c evaluated separately for the two sides is about 10 in 
both cases. The motion of the five atoms that displace the most is shown 
schematically relative to a reference atom (black), (b) The displacement 
of two (left) and three (right) atoms d as a function of the frame number 
j is shown schematically for a hypothetical pathway. The rearrangement 
on the left is more cooperative because two atoms move together over a 
longer region of the path. However, the current definition of N c does not 
distinguish between these two cases. 

6. N c as a function of N p calculated for the 57 most cooperative rearrangements 
from the LJ75 pathway database. 

7. N c as a function of N p calculated from four pathway databases for LJ 13 , 
LJ75 and BLJ256 systems. Due to the large number of data points we 
employ a density plot, where the darkest shading corresponds to the highest 
concentration of points. Outlying points are connected to illustrate the 
boundaries of the data area. The two BLJ databases are taken from Ref. 33 

8. Average barrier as a function of N p and N c calculated for the same LJ 13 , LJ75 
and BLJ 2 56 pathway databases used in Fig. The indices were calculated 
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separately for the two sides of each path. In this case the darkest shading 
corresponds to the highest barriers. 

9. N c as a function of N p calculated for the LJ75 pathway database. For each 
pathway we conducted DNEB calculations 8 assuming prior knowledge of 
the path. Every DNEB calculation employed 3s images and 90s iterations 
in each case, where s is the integrated path length. Out of 31,342 DNEB 
runs 25, 158 yielded a connected pathway while the rest did not (FAILED). 
Connected pathways are classified further as one-step pathways involving 
the correct transition state (OK) or an alternative transition state (ALT), 
or as multi-step pathways (MULTI), which involve more than one transition 
state. For each set of data points best fit straight lines obtained from linear 
regression are shown and labelled appropriately. 

10. The average value of N c for LJ 13 pathway databases as new paths are added. 
6-atom correlated perturbations (top curve) and random perturbations of 
every degree of freedom (bottom curve) were used to produce starting points 
for refinement by eigenvector-following. 17,20-22 ' 49 Average values were cal- 
culated every time 100 new pathways were added. 
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